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ABSTRACT 

A randomly fluctuating longitudinal pressure gradient of a non- Gaussian 
form and arbitrary power spectrum will cause a randomly fluctuating velocity 
to be superimposed on the steady incompressible flow in a channel. Pressure- 
gradient and velocity correlations, frequency response functions and system 
power loss are calculated. Numerical random pressure gradient signals 
were generated using Fourier series expansion, with random picking of dis- 
crete Fourier spectrum values and a nonlinear no memory transformation. 
Numerical values of the velocity signal were then calculated by linear trans- 
formation of the pres sure -gradient signal. Pressure -gradient and velocity 
signals were compared for difference in amplitude, frequency, time lag, and 
probability distribution functions. 

INTRODUCTION 

The study of random fluctuations in engineering system's is receiving 
increasing attention. This is because practically all physical systems have 
some randomness which at times can cause the system behavior to be dras- 
tically different from the steady case. 

The analysis of these stochastic problems is difficult, and usually only 
certain statistically averaged properties are obtained, such as the auto- 
correlation, which give only part of the information concerning the system . 
With the development of the high speed electronic computor it has become 
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feasible to study the random behavior of systems numerically, especially 
since the development of the fast Fourier transform algorithm which is a 
numerical method of rapidly taking Fourier transforms. This numerical 
calculation with random signals allows any statistical property of the signal 
to be calculated e.g. probability distribution function, by the appropriate 
numerical techniques. 

The present analysis describes the effect of a randomly fluctuating 
pressure gradient on the flow of an incompressible fluid in a channel (fig. 1) 
both analytically and numerically. A pressure -gradient signal having an 
arbitrarily specified probability distribution and power spectrum is gener- 
ated numerically. From this the velocity signal can be calculated. The 
velocity signal can then be analyzed for its probability distribution function, 
power spectrum and any other characteristic of interest. Analytical results 
are computed where possible and the numerical and analytical results are 
compared. 

In a previous analysis (ref. 1) this same flow problem was studied for 
the case where the pressure gradient fluctuations were restricted to a 
Gaussian (normal) probability distribution. Since the fluid velocity was 
linearly related to the pressure gradient, the velocity fluctuations would also 
be normally distributed. Therefore all the velocity characteristics would be 
given by the autocorrelation function of the velocity. However, for the non- 
Gaussian pressure gradient the velocity distribution is also non-Gaussian 
and cannot be characterized by the autocorrelation function alone. By use of 
a numerical method however the velocity signal is generated numerically and 
all its statistical characteristics can be found by the appropriate statistical 
analysis of the signal. 

Tien and Lienhard studied a similar problem for pipe flow (ref. 2) and 
included the effect of turbulent eddy diffusivity in Ref. 3. However only the 
autocorrelation results were given. This would not be sufficient to fully 
characterize a non-Gaussian system. 

As shown in Ref. 1 it is possible to generate numerically a Gaussian 
signal with a given power spectrum. Then by a nonlinear, no memory trans- 
formation this initial signal can be transformed into a pressure gradient sig- 
nal with the desired non-Gaussian probability distributions and associated 
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power spectrum. From this, through a linear transformation the velocity 
signal can be found and analyzed for its statistical properties. This process 
is outlined in Fig. 2. 


ANALYSIS 


In the present case we are Interested in generating a random pressure - 
gradient signal with a non-Gaussian probability distribution. To do so we 
must first generate a Gaussian signal then transform it to the desired signal. 
The method used is illustrated in Fig. 2. A Gaussian signal w(t) with a 
specified spectrum is subjected to a suitable nonlinear no memory transfor- 
mation such that a pres sure -graidient signal y is obtained having the desired 
probability density function p . This transformation can be written as 
y - g(w). We can relate the distribution functions of y and w by the well 
known relation 



The power spectrum of the generated signal S TT _„(w) must be specified 

Ww 

so that the nonlinear transformation, y = g(w), will have the desired power 

spectrum S . 

yy 

The autocorrelation of the generated w signal, which is the Fourier 
transform of the power spectrum, is obtained in terms of the desired pressure- 
gradient signal autocorrelation as follows. 


Generated Signal Autocorrelation 


We can write the pres sure -gradient signal autocorrelation as follows 
(ref. 4, pg. 51) 
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p is the normalized autocorrelation R /a„. and a. 
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is the variance 


It has been shown in Ref. 4 (p. 52) that the R can be written as an 


algebraic function of p Then for a known Gaussian signal, 
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transformation y = g(w); the autocorrelation p of the input Gaussian 

signal that, after the nonlinear transformation to the pressure signal, will have 

the desired autocorrelation R can be calculated by iteration. 
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A simpler method of calculating p is given in Ref. 5 for the case 
where the nonlinear transformation between 
specified form, 


y and w is given in the 
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The K and a are adjustable parameters that allow a wide variety of 
relationships between y and w to be approximated. 

The plot of y versus w given by Eq. (4) are plotted in Fig. 3. 

Although all values of w are possible the values of y are restricted between 
±0. 5/K. The value of K thus determines the cutoff amplitude of y. For 
a = 0 all values of w transform to ±0. 5/K giving a rectangular wave for y. 
As a becomes large, Eq. (4) reduces to Ky-*w/(n/27r ac^) so that yc cw. 

From Eqs. (1) and (4) we can write the distribution of y as 
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Eq. (5) can be evaluated by finding w for a specified value of y using 
Eq. (4). Then p /K can be evaluated from Eq. (6). The plot of p^/K 
versus yK is shown in Fig. 4. Significantly different probability distributions 
p can be obtained for different values of a and K. 

y 

If y = g(w) y given by Eq. (4) 5 is substituted into Eq. (2), then as shown 
in Ref. 5, R ; is given by 
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This can be evaluated at r = 0 to give 
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We can now write 
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which can be inverted to 


Pww = ( 1 + a ) sin 
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This gives the needed Gaussian autocorrelation p so that the r 
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transformation given by Eq. (4) will give the desired pressure gradient 
autocorrelation p . 
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Pressure Gradient Autocorrelation 


A common form of random fluctuation is called a stationary Gaussian 
Markoff Process. This gives an exponential form of the pressure gradient 
autocorrelation (ref. 6, pg. 215). 
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where the rate of fluctuation A and the time span T have been nondimen- 
sionalized by 
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where v is the fluid kinematic viscosity and d is the distance across the 
channel. The X is a dimensionless measure of the fluctuation rate. It can 
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also be considered the inverse of the dimensionless characteristic 
of the autocorrelation. The larger the X value, the greater the fluctuation 
rate of the pressure gradient signal. The fluctuation rate can be character- 
ized by the number of crossings of the o value of the signal per unit time. 

9 / 2 \ 

The quantity a is the mean square value of the pressure fluctuations \y f . 
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It also can be considered the integrated power spectrum of the pressure signal 
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where S / is the power spectrum and f is the frequency. The larger the ay, 

y y 

value, the larger the magnitude of the pressure fluctuations. 

The autocorrelation of the generated signal from Eq. (9) for the case of 
the stationary markoff pressure gradient signal becomes 


Pww = (1 + a) sin 
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The dimensionless autocorrelations and p , are shown in Fig 
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Since the power spectrum S is the Fourier transform of the auto - 
correlation, R^, we can write the dimensionless power spectrum of the pres- 
sure gradient as o 
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where co is the angular Fourier frequency, 27rf. The generated signal’s 
dimensionless power spectrum is given by 
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in a series give's, for the present pressure gradient 
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Generation of a Gaussian Random Signal by Model Sampling 

The Gaussian random signal w(t) can be obtained as shown in 
writing the signal as an inverse discrete Fourier Transform, 


rief. 1 by 
where n 
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are integer values. It follows that the discrete Fourier spectra can be con- 
sidered a Gaussian random variable with a zero mean and a variance that 
can be related to the power spectrum. We can use a random sampling pro- 
cedure to pick the discrete Fourier spectrum values from the appropriate 
distribution. Then by use of "fast Fourier transforms" we can obtain the 
signal values at discrete points. This detailed procedure is as follows. 

As shown in Ref. 1 the power spectrum, S , of a signal w which has 

a Fourier spectrum w, , can be written as 

x co r / > ”i 

(17) 

Where w* is the complex conjugate of the Fourier spectrum. The w a.R 
and w^j are the real and imaginary components of the Fourier spectrum. 

The term (. ) denotes ensemble averages. Thus the power spectrum can. be 
related to the average of the squares of the Fourier spectrum of the random 
signal. 

Discrete Fourier transform . - When a signal is to be analyzed on a 
digital computer it is the discrete Fourier spectrum rather than the con- 
tinuous spectrum that must be considered. If the Fourier spectrum w 
is band limited so that w^ ~ 0 for f > f ^ we can write (ref. 7, pg. 152) 
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Thus knowing w at discrete values erf t allows w w to be evaluated from 

Eqs. (18) and (20). Then by taking the inverse continuous Fourier transform 

of w we obtain w at all times. Since we assumed t = n/f = nAt we can 
co p 

see that 



cok \ e i2?rkn/N 


k^/2 ' P / 

So that Eqs. (23) and (24) form a discrete Fourier transform pair that are 
approximations to the continuous Fourier transform pair. 

Fast Fourier transforms . - A numerical method has been developed, 
called the fast Fourier transform (ref. 8), for evaluating discrete Fourier 
transforms in a very efficient manner on a digital computer. One imple- 
mentation readily available is the IBM -scientific subroutine package caiiec 
HARM /D HARM , Many others are available. 


Model sampling to generate Gaussian signal . - We can expand w^ k 
in Eq. (24) into real and imaginary parts. Since w n is a real function 
the imaginary term must be zero and so w must be an even function 
around n = 0 while w^j is an odd function around n = 0. Then we can write 
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It has been shown in Ref. 1 that for a Gaussian random process w wn ^ 
and w j are normally distributed, independent, random variables with 
zero mean and variance given by 
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n = 0 


We can generate the values of w ^ and w wn j by first writing their 
joint Gaussian probability distribution 
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Following the procedure in Ref. 1 we can randomly sample values from t! 
above distribution with the equations 
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where R 
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and R f are two different random numbers picked from a uniform 


distribution that is constant between 0 and 1. These values can readily be 
generated on a digital computer. Recalling that w wnR is an even function 
and w^ n - is an odd function, in the particular fast Fourier program used 
(HARM/DHARM), we need only find w wn ^ and w^j from n = 0 to N/2. 

Then w^ ( N / 2 ) + -l, R = w co, (N/2)-l, R» etc ‘ Simllarly w co, (N/2)+l, I = 

-W co (N/2)-l I’ etc ‘ The value of w co 0 I is e T ual to zero. 
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After finding the values of we can substitute them into Eq. (2! 

and using the fast Fourier transform routine find the values of w, . 

The sample mean w and the sample variance are then found ir< 


the w k values 


K 


With this the value of w, can be narmalizec 

w)/(wV /! . 


Then the distribution and autocorrelation are fount 


w k " xv# k 

by using available subroutines on the values of w^. 

The example case is carried out using the power spectrum S^,. given 
by Eq„ (16). The normalized signal w 7 is plotted in Fig. 6 for a X value 
of 1. 


In Fig. 7 the cumulative distribution function of w 7 obtained numeri- 
cally is plotted along with the Gaussian cumulative distribution function with 
good agreement. This shows that the w 7 is normally distributed as expected. 

The numerical results for the autocorrelation of the w 7 signal is shown 
in Fig. 5 and compares well with its theoretical result, although the agree- 
ment seems to fall off for large time lag t. 

Generation of y signal from w signal . - To transform the w signal 
into the pressure gradient signal y we use Eq. (4) which allows us for each 
w R (nAt) to calculate a y n (nAt) value for each n from 0 to N. This gives 
the pressure gradient signal with the desired probability distribution and 
autocorrelation. This signal is normalized to y 7 and shown in Fig, 6. 

Since the normalized pressure gradient signal cannot have any values larger 
than (0. 5/a^,) we can see this cutoff in Fig. 6. It can be seen from Fig. 4 
that the pressure cutoff is at lower values of y 7 for smaller values of X 
since the normalizing factor cr increases for decreasing a, see Fig. 4. The 
cumulative distribution of the y ? signal is shown in Fig. 7 where it can he 
seen to agree with the analytical result. The autocorrelation of the y 7 signal 
is shown in Fig. 5 and also agrees with the analytical result. 


Relation Between Velocity Signal and Pressure Signal 

This procedure follows the analysis in Ref. 1. To calculate the fluid 
velocity signal from the pressure gradient signal we must first solve the 
momentum equation. For viscous incompressible flow between parallel plates 
with constant properties in the fully developed flow region (see fig, 1) the 
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momentum equation can be written as 
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We can let the velocity and the pressure gradient consist of a steady part and 
a nonsteady part U = U_ + U, ; P = P + P, . Normalizing as follows 
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The solution for v as shown in Ref. 1 is given by 
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where M = m77. The Fourier transform of Eq. (32) is given by 
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and Y m = 4 sin My/M. The H is sometimes called the system transfer func- 
tion. We can write the real and imaginary parts of H as 
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The gain factor of the system function is given by Jh| = (H-^ + Hj) '" 7 and the 
phase factor of the system function by c p = tan - ^ (Hj/Hp). The gain and phase 
factors are plotted in Fig. 8. As can be seen in the figure, for small values 
of angular frequency c o the system frequency transfer function approaches the 
distortionless limit |h| = | h]^; cp^ = - cot ^ to give 
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Then taking inverse Fourier transforms of Eq. (34) 
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This shows for very low frequency, co — 0, the velocity signal duplicates the 
pressure signal but is a fraction | H J ^ of the pressure signal amplitude. The 
velocity signal lags the pressure signal by a time increment t r We can find 
the limiting value of t^ as co — 0 as 0.0781. 

The gain factor decreases in magnitude close to the wall. This shows 
that the magnitude of the velocity fluctuations close to the wall is smaller than 
near the center of the channel. The percent decrease in the gain factor at the 
higher frequencies is larger near the center of the channel so that there will 
be relatively more damping of the higher frequency velocity fluctuations near 
the center of the channel than near the wall. The phase factors are smaller 
near the channel wall than near the center showing that there is a smaller 
signal lag of the velocity signal behind the pressure signal closer to the wall. 

The cross correlation of the pressure and velocity fluctuations can be 
obtained from Eq. (32) as 

/% oo ^ OO 

R yv (r) = (y(t)v(t + t )) =J h(0, y) (y(t)y(t + t - 0))d 6 =1 h(6y)R yy (T - 8)d6 

J 0 J 0 ( 38 ) 
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where R is the time lag. Since the cross spectrum is the Fourier trans- 
form of the cross correlation, by using the convolution theorem we can write 

S yv («) - (39) 

We can find the velocity autocorrelation similarly as 
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Then we can write the power spectrum as 
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The pressure -velocity cross correlation for the present case can be 
written using Eqs. (10), (33), and (38) as 
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In a similar manner we can calculate the velocity autocorrelatior 
from Eqs. (10), (33), and (40) as 


OO OO 


R (r) = cT 

W v ' y 


m=l, 3 m=l, 3 


y m y m- 


(M’ 2 - X) 


— r 

) i 


e~ Xr /(\+M‘ 


- ^2X/[x+M , 2 1) e" Mt T /(M 2 +M ,2 )J (43) 

The normalized power spectrums are plotted in Fig. 9. The pressure gradient 
spectrum from Eq. (14) for the stationary markoff process can be seen to have 
a lower frequency cutoff for lower values of X or slower rates of fluctuations. 
The normalized velocity spectrums given by Eq. (41) are also shown and for 
the lower values of X are close to the normalized pressure gradient spectrums 
and not strongly dependent on distance across the channel. However at the 
higher fluctuation rate, or large X, the power spectrum of the velocity sig- 
nal have more of the power concentrated at the lower frequency then the power 
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spectrum of the pressure gradient signal. 

The mean square velocity, (v 2 ), can be obtained from Eq. (43) since 

f» 4-00 » -}-oo 

(v 2 ) = R vy (o, y) = / S yv df = I | H(o>) 1 2 S yy (a))df (44) 

J- oo J - oo 

The values of ^v 2 ^ across the channel for different values of X are shown 
in Fig. 10. 

Since the gain factor is smaller near the wall (y = 0. 1) compared to that 
at the centerline (y = 0.5), the mean square velocity is smaller near the wall 
as compared to the centerline. The mean square velocity goes to zero at the 
wall (fig. 10) as would be expected from the zero wall velocity boundary 
condition. The gain factor is greatest at the lower frequency and has a high- 
frequency cutoff (cf. fig. 8). Therefore for a given mean square pressure 
gradient value the maximum value of ^v 2 j is found when the pressure 
gradient spectrum is a delta function at small w. The minimum value of 
(v 2 ) occurs if the pressure spectrum is a delta function at large value 
of w. 


Power Dissipation Due to Velocity Fluctuations 


It is desirable to know how much additional power must be supplied to the 
flow system to maintain the same average flow rate for the fluctuating flow as 
compared to the steady flow case. As shown in Ref. 1, the external rate of 
work done by the pressure force w g is equal to the sum of the rate of increase 
in kinetic energy K g and the rate of dissipation of energy due to internal 
friction w^. If we then take the ensemble average, the rate of change of the 
kinetic energy term will be zero because the process is stationary. We then 
have the external rate of work done by the pressure force w g equal to the 
rate of dissipation of energy due to external friction w^ 


where 
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U is the integrated value across the channel, 



( 45 ) 
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To find the rate of dissipation of energy due to internal friction we can 
then write 

w f = - (dp) = - T3 S P S - (u t P t ) - * (f „ - (u t P t ) 

If we let w, be the rate of dissipation due to steady flow, the ratio of 
I y s 

increased dissipation due to the unsteady flow over the steady flow dissipa- 
tion can be written as 
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Integrating Eq. (42) over y and substituting the result Eq. (51) give 

OO 

(Wf - w^ g )/Wj s = 96cr^ jjT 1/M^(M^ + X) (48) 

m=l, 3, 5 

This result is shown in Fig. 11 as a function of the fluctuating rate param- 
eter X. As X becomes very large, that is, as fluctuations of the pressure 
gradient become very rapid, the frictional power loss reduces to the steady 
power loss. The frictional power loss increases with smaller values of X, 
that is, with slower fluctuations of the pressure gradient. For a given mean 
square pressure gradient the maximum rate of power dissipation that can be 
achieved can be seen from Eq. (47) to be when is a delta function at 
the maximum value of H which is at small oo. The minimum power dissipa- 
tion occurs when S is a delta function at the minimum value of E: i. e. , 

yy 

oo large. 


Calculating Velocity Signal From Pressure Signal 


We can write the velocity signal as before (eq. (24)) 
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From Eq, (34) we can write 

v con 
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The values of y^ n can be obtained by taking the fast Fourier transform 
of y^ which were calculated previously. Then by use of Eq. (50) we can 
obtain the values for v un » The values of v n are then obtained using the 
inverse fast Fourier transform. From these values the normalized values 
v^ can be found. These are plotted in Fig. 12. When compared to the pres- 
sure signal it can be seen there is rounding of the sharp cutoff of the pressure 
signal. Fig. 7 shows the cumulative distribution of the velocity signal 
which can be seen to be similar in shape to the cumulative distribution of the 
pressure signal but tending towards the Gaussian cumulative distribution. The 
autocorrelation of the velocity signal is shown in Fig. 13 to agree fairly well 
with the theoretical results given by Eq. (43) but diverges somewhat for large 
lag times t. 


SUMMARY AND CONCLUSIONS 

The effect of imposing a randomly fluctuating pressure gradient with a 
non-Gaussian probability distribution and an arbitrary spectrum on si fluid in 
a channel causes a randomly fluctuating velocity component to be superimposed 
over the steady flow. The mean value of the velocity fluctuation is zero. The 
mean square value of the velocity fluctuations which is a measure of the am- 
plitude of the velocity fluctuations is highest in the center of the channel and 
reduce to zero at the wall. Also the slower the rate of the pressure gradient 
fluctuations the greater the amplitude of the velocity fluctuations. 

The pumping power loss was increased by the flow fluctuations;, the slower 
pressure gradient fluctuation rate giving the larger power loss. 

The non-Gaussian random pressure gradient signal with a stationary mark- 
off power spectrum was obtained numerically by first generating a Gaussian 
signal with a specified power spectrum. This was generated by first choosing 



17 


random values of the coefficients of the discrete Fourier spectrum repre- 
sentation from a Gaussian distribution whose variance was obtained from the 
specified power spectrum. Taking fast Fourier transforms of the discrete 
Fourier spectrum gave discrete values of the Gaussian signal. This signal 
was then transformed through a nonlinear, no memory, relation to the de- 
sired non-Gaussian fluctuating pressure gradient signal. The velocity signal 
was obtained by first taking the fast Fourier transform of the pressure gra- 
dient signal and multiplying the resulting discrete Fourier spectrum by the 
system frequency response function. This gave the discrete velocity Fourier 
spectrum which was then fast Fourier transformed to the velocity signal. The 
numerical signals were then analyzed to obtain their correlations and frequency 
distribution functions. 

The normalized velocity signals were close in value to the normalized 
pressure gradient signal but had a -dimensionless time lag vT/d of about 
0.0781 for the lower pressure gradient fluctuation rates (small X), At the 
higher rates of fluctuation of the pressure gradient signal (large X) less of 
the high frequency fluctuations pressure gradient signal appear in the veloc- 
ity signal. More so near the center of the channel rather than near the wall. 
Also the time lag of the velocity signal behind the pressure gradient signal is 
larger near the center of the channel then near the wall. The velocity signal 
tended to overshoot the upper limit of the pressure gradient signal. 

The nondimensional cumulative probability distribution for the velocity 
signal is very close to the cumulative probability distribution for the pressure 
gradient at the low fluctuation rate (small X). However at the higher values 
of X the normalized cumulative probability distribution for the velocity 
approached the Gaussian cumulative probability distribution. 
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APPENDIX - SYMBOLS 

complex coefficients 
spacing between parallel plates 
expected value 
frequency 

frequency increment, 1/T 
maximum frequency 

Fourier transform, inverse Fourier transform 
frequency response function 
gain factor 

weighting function or impulse response function 
unit imaginary number 
m-ir; m®ir 

number of samples 

pressure gradient 

probability distribution function 

correlation (x, y) 

spectrum 

time 

period of record 

2 

dimensionless time, vT/d 
time between sample 
fluid velocity 

2 

steady dimensionless velocity, -vplj /d P 

s s 
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2 

v transient dimensionless velocity, -I'pUt/d P 

s 

w generated Gaussian signal 

Y coordinate across channel 

Y m 4 sin My/M 

y dimensionless coordinate across channel, Y/d 

a parameter (eq» (4)) 

p constant; sin” (1/1 + a ) 

y dimensionless pressure gradient, P./P 

L S 

0 parameter given by Eq. (7) 

A measure of rate of fluctuation of pressure signal 

A dimensionless measure of rate of fluctuation of pressure 

signal, M 2 /v 

v kinematic viscosity 

p_„ normalized autocorrelation R _ /R (0) 

2 

a variance 

r dimensionless time difference, tg - t^ 

cp phase factor 

normalized power spectrum, S/R(0) 
co angular frequency, 2 77 f 

( ) sample average 

( ) ensemble average 

* complex conjugate 

’ normalized variable 
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Subscripts: 

F fundamental frequency or lowest frequency, 1/T 

H related to frequency response function 

I imaginary 

R real 

s steady component 

t fluctuating component 

y related to pressure signal 

oo Fourier spectrum 
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Figure 4. - Probability distribution of pressure gradient 
signal. 
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Figure 6. - Normalized generated Gaussian signal and pressure gradient 
signal. 
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Figure 7. - Normalized cumulative distributions. 
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Figure 8. - System transfer function. 
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Figure 10. - Normalized mean square fluctuating velocity com- 
ponent across channel. 



Figure 11.- Percentage increase in power dissipation ratio due to 
flow fluctuations. 
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Figure 12. - Normalized pressure gradient and velocity signal. 



;-6393 


ANALYTICAL CORRELATION NUMERICAL 

PRESSURE GRADIENT p™ a 

vaOCITY Q, a- 1; □, a= 0. 5 



DIMENSIONLESS TIME LAG, Xt 


Figure 13. - Dimensionless autocorrelation of pressure gradient and velocity 
signal. 
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